Computing intersections with polygons¶

Binder IPYNB HTML

Clipping and intersection functions can be used to extract trajectory segments that are located within an area of interest polygon.

In [1]:
import pandas as pd
import geopandas as gpd
from geopandas import GeoDataFrame, read_file
from shapely.geometry import Point, LineString, Polygon
from datetime import datetime, timedelta
import movingpandas as mpd

import warnings
warnings.filterwarnings('ignore')

mpd.show_versions()

MovingPandas 0.9.rc2

SYSTEM INFO
-----------
python     : 3.7.12 | packaged by conda-forge | (default, Oct 26 2021, 05:37:49) [MSC v.1916 64 bit (AMD64)]
executable : E:\Anaconda\envs\mpd-ex\python.exe
machine    : Windows-10-10.0.19041-SP0

GEOS, GDAL, PROJ INFO
---------------------
GEOS       : None
GEOS lib   : None
GDAL       : 3.2.1
GDAL data dir: None
PROJ       : 7.2.0
PROJ data dir: E:\Anaconda\envs\mpd-ex\Library\share\proj

PYTHON DEPENDENCIES
-------------------
geopandas  : 0.10.2
pandas     : 1.3.5
fiona      : 1.8.18
numpy      : 1.21.5
shapely    : 1.7.1
rtree      : 0.9.7
pyproj     : 3.1.0
matplotlib : 3.5.1
mapclassify: 2.4.3
geopy      : 2.2.0
holoviews  : 1.14.6
hvplot     : 0.7.3
geoviews   : 1.9.2
In [2]:
gdf = read_file('../data/geolife_small.gpkg')
traj_collection = mpd.TrajectoryCollection(gdf, 'trajectory_id', t='t')

Clipping a Trajectory¶

In [3]:
help(mpd.Trajectory.clip)

Help on function clip in module movingpandas.trajectory:

clip(self, polygon, point_based=False)
    Return trajectory segments clipped by the given polygon.
    
    By default, the trajectory's line representation is clipped by the
    polygon. If pointbased=True, the trajectory's point representation is
    used instead, leading to shorter segments.
    
    Parameters
    ----------
    polygon : shapely Polygon
        Polygon to clip with
    point_based : bool
        Clipping method
    
    Returns
    -------
    TrajectoryCollection
        Clipped trajectory segments
In [4]:
xmin, xmax, ymin, ymax = 116.365035,116.3702945,39.904675,39.907728
polygon = Polygon([(xmin,ymin), (xmin,ymax), (xmax,ymax), (xmax,ymin), (xmin,ymin)])
polygon_gdf = GeoDataFrame(pd.DataFrame([{'geometry':polygon, 'id':1}]), crs=31256)

my_traj = traj_collection.trajectories[2]
intersections = my_traj.clip(polygon)
    
print("Found {} intersections".format(len(intersections)))

Found 1 intersections
In [5]:
ax = my_traj.plot()
polygon_gdf.plot(ax=ax, color='lightgray')
intersections.plot(ax=ax, color='red', linewidth=5)
Out[5]:
<AxesSubplot:>

Clipping a TrajectoryCollection¶

Alternatively, using TrajectoryCollection:

In [6]:
clipped = traj_collection.clip(polygon)
clipped
Out[6]:
TrajectoryCollection with 2 trajectories
In [7]:
clipped.plot()
Out[7]:
<AxesSubplot:>

Computing intersections¶

In [8]:
help(mpd.Trajectory.intersection)

Help on function intersection in module movingpandas.trajectory:

intersection(self, feature, point_based=False)
    Return the trajectory segments that intersects the given feature.
    
    Feature attributes are appended to the trajectory's DataFrame.
    
    By default, the trajectory's line representation is clipped by the
    polygon. If pointbased=True, the trajectory's point representation is
    used instead, leading to shorter segments.
    
    Parameters
    ----------
    feature : shapely Feature
        Feature to intersect with
    point_based : bool
        Clipping method
    
    Returns
    -------
    TrajectoryCollection
        Segments intersecting with the feature
In [9]:
polygon_feature = {
    "geometry": polygon,
    "properties": {'field1': 'abc'}
}
In [10]:
my_traj = traj_collection.trajectories[2]
intersections = my_traj.intersection(polygon_feature)
intersections
Out[10]:
TrajectoryCollection with 1 trajectories
In [11]:
intersections.to_point_gdf()
Out[11]:
id sequence trajectory_id tracker geometry intersecting_field1
t
2009-02-04 10:43:04.222205 3225 773 3 2 POINT (116.36799 39.90467) abc
2009-02-04 10:43:05.000000 3226 774 3 2 POINT (116.36798 39.90471) abc
2009-02-04 10:43:07.000000 3227 775 3 2 POINT (116.36795 39.90480) abc
2009-02-04 10:43:09.000000 3228 776 3 2 POINT (116.36793 39.90487) abc
2009-02-04 10:43:11.000000 3229 777 3 2 POINT (116.36794 39.90495) abc
... ... ... ... ... ... ...
2009-02-04 10:48:14.000000 3390 938 3 2 POINT (116.36535 39.90589) abc
2009-02-04 10:48:15.000000 3391 939 3 2 POINT (116.36525 39.90589) abc
2009-02-04 10:48:16.000000 3392 940 3 2 POINT (116.36516 39.90589) abc
2009-02-04 10:48:17.000000 3393 941 3 2 POINT (116.36507 39.90589) abc
2009-02-04 10:48:17.399995 3393 941 3 2 POINT (116.36504 39.90589) abc

170 rows × 6 columns

In [ ]: